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Abstract 

This paper develops a structure-preserving numerical integration scheme for a class 
of higher-order mechanical systems. The dynamics of these systems are governed by 
invariant variational principles defined on higher-order tangent bundles of Lie groups. 
The variational principles admit Lagrangians that depend on acceleration, for exam- 
ple. The symmetry reduction method used in the Hamilton-Pontryagin approach for 
developing variational integrators of first-order mechanics is extended here to higher 
order. The paper discusses the general approach and then focuses on the primary ex- 
ample of Riemannian cubics. Higher-order variational integrators are developed both 
for the discrete-time integration of the initial value problem and for a particular type of 
trajectory-planning problem. The solution of the discrete trajectory-planning problem 
for higher-order interpolation among points on the sphere illustrates the approach. 
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1 Introduction 

After introducing the general background for higher-order mechanics and variational integra- 
tors, we state our goals and summarize the content of the paper. 

1.1 General background 

The problem treated in this paper fits into a classical type of problem in control theory called 
trajectory planning, or interpolation by variational curves. The task in this type of problem 
is to find an optimal curve that interpolates through a given set of points (or configurations) 
lying in a manifold. The configuration manifold is specific to the application. In this paper, we 
shall study the example of trajectory planning for tracking a rigid body through a prescribed 
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sequence of orientations. The control of a rigid body summons either the group SO (3) of 
rotations in IR 3 , or the semidirect-product group SE(3) = SO (3) X M. 3 of three-dimensional 
rotations and translations in Euclidean space. Such trajectory-planning problems are relevant 
in numerous applications, for example, in aeronautics, robotics, computer-aided design, air 
traffic control, biomechanics and more recently, computational anatomy. 

Higher order mechanics. Some trajectory-planning applications may require the optimal 
trajectories to possess a certain degree of smoothness. This requirement summons variational 
principles that depend on higher-order derivatives of the interpolation path, such as acceler- 
ation (rate of change of velocity) or jerk (rate of change of acceleration) etc. Properties of 
these higher-order variational principles, including their Hamiltonian formulation and their 
symmetry reduction, have been studied in [dLR85, BdD05, GBHR11, GBHM+11, CdDll]. 
An interesting example of such a higher-order variational principle was introduced in [GK85] 
and [NHP89]. This example requires the minimization of the mean-square covariant ac- 
celeration and leads to curves called Riemannian cubic splines. These curves generalize 
the familiar cubic splines of Euclidean space to the context of Riemannian manifolds. The 
mathematical theory of Riemannian cubics and their higher-order generalizations has been 
developed in a series of papers [CS95, CSC95, CSC01, Noa04, Noa06a, GGP02]. Applica- 
tions to computer graphics, spacecraft control and computational anatomy are discussed in 
[PR97, ZKC98, HB04, VT11], amongst others. We refer to [Pop07, MSK10] and [Noa06b] for 
extensive references and historical discussions concerning Riemannian cubics, their higher- 
order generalizations, and related higher-order interpolation methods. 

Many properties of first-order mechanics remain relevant in the higher-order setting. In 
particular, curves that solve Hamilton's principle for a given non-degenerate Lagrangian in 
first-order mechanics have various interesting geometric features. For example, the corre- 
sponding Hamiltonian and symplectic structure, both defined on the cotangent bundle of the 
configuration manifold, will be conserved along these solutions. In addition, if the Lagrangian 
is invariant under the action of a Lie group then the associated momentum map is also con- 
served (we refer to [MR03] for a detailed discussion). In fact, much of this structure has been 
extended to higher order [dLR85, GBHM + 11, GBHR11, CdDll]. Moreover, in first order 
mechanics much has also been achieved in bringing these properties to the numerical side, 
where time is discretised and one aims to approximately solve the variational equations. The 
numerical schemes oriented towards preserving such properties form the topic of geometric 
integration, which is a relatively new subject for higher-order mechanics [HLW06, CJdDll]. 

Variational integrators. In designing a geometric integrator, one must face the fact that 
all three properties (momentum preservation, symplecticity and energy preservation) cannot 
be achieved simultaneously by methods with a fixed time step [KM099, Mar92, HLW06]. Our 
focus is on the symplectic and momentum preserving schemes that follow from discretising 
Hamilton's variational principle. The resulting variational integrators are well established 
| Cad 70, Mae82, WM97] and the discrete analogues of many of the familiar properties implied 
by Hamilton's principle, such as Noether's theorem, have been established [Log73]. This 
progress has made variational integration a convenient approach for numerics associated 
with mechanics [Lee87]. Despite lacking energy preservation, variational integrators that are 
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symplectic benefit from the favourable energy behaviour and long time reliability properties 
afforded by symplectic methods for Hamiltonian systems [Row91, BG94, Rei99, HLW06]. The 
background and breadth of variational integration for Lagrangian mechanics on manifolds is 
reviewed by Marsden and West [MW01]. 

For our purposes here, the preservation of Lie group structure is also relevant. This is- 
sue has been addressed for many types of geometric integration [LS94, HLW06, IMKNZOO, 
AKW93]. Variational integrators for first-order mechanics on Lie groups and the correspond- 
ing discrete Euler-Poincare and Lie-Poisson equations have been discussed, for example, in 
[MPS99, BS99, LLM07]. The prevalence of Lie groups in applications means these integrators 
have been implemented in various fields, including computer graphics, integrable systems, 
solid mechanics and quantum mechanics [KYT + 06, MV91, OS99, JN97]. 

An interesting variant in the approach to these Lie group integrators was introduced in 
[BRM09], where discretisation was applied to the first-order Hamilton-Pontryagin (HP) vari- 
ational principle [H0IO8, YM06], rather than to the more familiar Hamilton's principle. The 
HP principle is of particular significance to the present investigation because it explicitly 
contains the definitions of the variables of interest for mechanics, position, velocity and mo- 
mentum, which are all elements of the Pontryagin bundle [BRM09]. This explicit dependence 
in the HP principle allows one to take unconstrained variations over the Pontryagin bundle. 
This property carries over to the discrete HP principle and simplifies certain aspects of dis- 
crete mechanics. For example, it removes the need to identify discrete Legendre transforms 
when discussing symplectic structure. This property greatly facilitates the move to discrete 
higher-order mechanics. 

1.2 Motivation 

Our primary motivation is to develop symplectic, momentum-preserving integrators for higher- 
order mechanics on Lie groups. The methods used to develop Hamilton-Pontryagin integra- 
tors in [BRM09] are particularly amenable to extension to higher order due to the explicit 
nature of the HP variational principle. We investigate the geometric properties of the result- 
ing numerical algorithms both for the initial value problem and for applications in trajectory 
planning. 

1.3 Main content of the paper 

The plan of the paper is as follows. 

After introducing the basic definitions that we will need, Section 2 develops a theorem 
that characterises Hamilton's principle for second order mechanics on Lie groups (anal- 
ogous to the first order principle) for Lagrangians with symmetry. We show equivalence 
with two higher-order versions of variational principles from first-order mechanics: (i) 
the reduced Hamilton's principle and (ii) the reduced Hamilton-Pontryagin (HP) prin- 
ciple [BRM09]. In this paper, our main focus is on the latter, in which Lagrange 
multipliers are used to enforce particular kinematic constraints. We use this princi- 
ple to derive some previously known results for continuous-time systems, including the 
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Fig. 1.1: Illustration of the trajectory-planning problem on the sphere discussed in Section 4. In this problem, 
one seeks a curve that passes near the target points (black dots) at prescribed times. The trajectory shown 
is generated from a path of rotations acting on Euclidean space. The variational description requires the 
minimization of a functional that measures both the amount of acceleration (measured in the group of 
rotations) along the curve and the amount of mismatch between the targets and the curve at the prescribed 
times. The discussion in Section 4 shows that the geometric properties of optimal curves are inherited by 
the discrete higher-order Hamilton-Pontryagin description of the problem. 



NHP equation of [NHP89] for Riemannian cubics. Continuing in Sections 2.3 and 2.4 
we describe the Hamiltonian theory and symmetry properties underlying higher-order 
variational problems on Lie groups. In particular, we will discuss how the flow map 
induced by the HP principle is symplectic and preserves momentum. 

In Section 3 we discuss discrete mechanics and obtain the discrete geometric counter- 
parts of the HP differential equations using a method similar to that in [BRM09]. The 
kinematic constraints are discretised using Runge-Kutta methods and their Lie group 
incarnations, the Runge-Kutta-Munthe-Kaas methods [IMKNZ ]. These are then 
incorporated by using Lagrange multipliers into a discrete version of the higher-order 
HP principle. We will find that the flow maps of the resulting algorithms inherit the 
symplectic momentum-preserving properties of the continuous-time systems in Section 
2. Having established the theory, we proceed in Section 3.2.2 with an implementation 
of it in the particular case of Riemannian cubics on 5*0(3) with a bi-invariant metric. 

Section 4 applies these discrete methods to a second order trajectory-planning problem 
in which the interpolating curve evolves by means of a group action. We first formulate 
the problem in the framework of the Hamilton-Pontryagin principle. As in previous 
work [GBHM + 11], we find that the momentum experiences kicks related in a simple 
fashion to the mismatch between the interpolating curve and the target points - a 
matter which can be presented particularly clearly in the HP framework. We show how 
the discretisation of the problem respects this geometric feature of solution curves. The 
section proceeds with a discussion of a numerical algorithm that allows us to calculate 
curves such as the one shown in Figure 1.1 for the trajectory-planning problem on the 
sphere. 
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In Section 5 we summarise and present a series of outstanding directions for further 
research. 

2 Continuous-time mechanics 

Here we develop the theory of higher-order mechanics in the formalism of the Hamilton- 
Pontryagin variational principle. After giving some general definitions, we introduce the 
fundamental dynamical equations. We then proceed to discuss Hamiltonian structure and 
momentum maps. 

2.1 Definitions 

We start by introducing the main mathematical objects used in the paper. Let G be a Lie 
group and k a positive integer, k > 0. The k th -order tangent bundle Tq : T^G — > G is 
defined as a set of equivalence classes of curves, as follows: Two curves ji : 1 1— > g%(t), i = 1, 2, 
are equivalent^ if and only if their time derivatives at t = up to order k coincide in any local 
chart. That is, g[ (0) = g% (0), for < I < k. The equivalence class of a curve 7 : t 1— > g(t) 
is denoted by [7]^)' or f° rma lly as (#(0), <7(0), . . .g^(0)). The set of all equivalence classes 
of curves based at g is written as Tg^'G. 

The k th -order Pontryagin bundle P^G is defined as [CdDll] 

P (k) G = T {k) Q XT(fc _ 1)G t*{T^G) 

Note that T^G = G, T (1) G = TG and P«G = TG x T*G. The projection map onto the 
second factor is denoted by pr2 : P^ 2 'G — > T*(TG). For any fixed real number a > and 
smooth manifold M, we write C(M) to denote the set of smooth curves 7 : [0,a] —> M, 
t^^t). 

The group structure of G affords the following trivialisations which make extensive use 
of the left multiplication maps, L g : G — > G, h i-> #/t. We have TG = G x q via f 9 1— > 
(g,TL g -iv g ) = (g,g- l v g ), and T*G = G x q* via a 9 1— >■ (g, (TL g )*a g ) = (g,g*a g ). Here 
we introduced the notation g~ 1 v g := TL g -iV g and := (TL g )*a g . In order to define 

the trivialisation of T^G let g(t) be a representative of [7]^) anc ^ ^(0 = 9~ X {^)9{^) ^ e the 
left-trivialised velocity. The trivialisation map T^G = G x k$ is given by 

WSJ,-* (»(o),e(o),e(o) ) ... ) e ( *- 1) (o))- 

These maps further induce the identification T* {T^ 1 ^ G) = G x [k — l)g x k$*. Accordingly, 
p(k)Q ^ Q x k$ x kQ*. Note that all of these trivialisations make use of the group multi- 
plication from the left. Right-trivialisations can be defined correspondingly using the right 
multiplication. 

From here on we shall consider second order mechanics, k = 2, and for later reference we 
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set the notation 

( 9l ( 1 i 1 )eGx2 6 -T( 2 >G 1 
(g,^fx,v)eGxQx2 Q *=T*(TG), 
(g,C,u,ii, i/)GGx2 x2 S * = P (2) G. 
We also define the projection II : P^G — > G onto the first factor. 

2.2 Variational principles 

A second-order Lagrangian L : T^G — > K. is called left-invariant if its left trivialisation 
I : G x 2g — > K. does not depend on the first entry. One then introduces the reduced 
Lagrangian £ : 2q — > R. The following theorem, the cornerstone of our approach, can easily 
be generalised to higher order, k > 2. 

Theorem 2.1. Let L : T^G — > R be a left-invariant Lagrangian with reduced Lagrangian 
£ : 2q — > K. and let g be in C{G). The following are equivalent: 

1. Hamilton's principle. The curve g{t) satisfies 5S(g) = for S : C{G) — > WL, 

S = [ L(g,g,g)dt, (2.1) 
Jo 

with respect to arbitrary variations 5g G T g C(G) which fixed g(0) , g(a), g(0), g(a). 

2. Reduced Hamilton's principle. The curve v(t) := g~ l {t)g{t) satisfies 5s(v) = for 
s : C(fl) ->• K, 

r 

s — £{v, v) dt, 
Jo 

with respect to variations of the form 5v = t) + &d v i], 77 G C(q) with 77(0) = 77(a) = 
77(0) = 77(a) = 0. 

3. Reduced higher- order Hamilton-Pontryagin (HOHP) principle. There exists a curve 7 
in P^G with projection II o 7 = g that satisfies 05^(7) = for Sh p ■ C(P^G) — > M, 

s hp = J e& u) + </i, g- l g + uj dt, (2.2) 

with respect to variations £7 G T 7 C(P^) with 5g(0) = 5g(a) = and <5£(0) = <5£(a) = 0. 

Equivalence of these variational principles follows from comparing the equations they 
imply. For illustration, we derive the equations governing solutions to the HOHP principle. 
For 7 G C(P^G) as above denote by 27 G T 1 C{P^ > G) an arbitrary vector over 7. That is, 
Zry : t (->■ (5g(t), 5£(t), 5u(t), 5/i(t), 5u{t)). We define 77 = g~ l 5g and, upon recalling Si^g^g) = 
rj + ad 9 -ig 77, we compute 



os/^ = ds fcp (2 7 ) = [(//, 77) + (MOlo + y o < ad g-i9/ i -A^> + (^-/ i - z> > ^ 
+ ( ^ ~ V ' 6U ) + 9 ^ ~~ + i-uj dt. 



(2.3) 
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The end-point terms will be of interest later, but for now we employ the HOHP principle 
with f](0) = r](a) = 5£(0) = 5£(a) = 0. We obtain the HOHP equations of motion, 

fi = ad^ /i, (2.4) 

81 . 5£ 
/,= --,,,= -, (2.5) 

9 = 9t, Z = u, (2.6) 

which agree with the ones derived, in [GBHM + 11], from the reduced Hamilton's principle. 
Equation (2.4) is usually called the second order Euler-Poincare equations. Equations (2.5) 
define the Ostrogradsky momenta fi and v. 



Riemannian cubics. The paper focuses on Riemannian cubics. These variational curves 
were introduced in [GK85] and [NHP89] in the general context of Riemannian manifolds. 
Here we consider the special case of cubics on Lie groups. 

Let G be a Lie group with a Riemannian metric d and associated norm functions ||.|| : 
T g G — > E. We define the isomorphism b : $j — > Q* by = d(r),£), for all (eg. Denote 

its inverse by Jj = b _1 . Consider the Lagrangian L : T^G — > R, 

(2.7) 



L(g,g,g) = \ 



D 

Dt 9 
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where D/Dt is the covariant derivative along curves, with respect to the Levi-Civita con- 
nection for the metric d. If the metric is left-invariant then L is left-invariant with reduced 
Lagrangian £ : 2q — > K, 

£(e,«) = i||«-ad € e t ||^ (2-8) 

where ad^ is the metric adjoint, that is, 

ad^r/ := (ad^r/ b )« , (2.9) 

for all £, r] e Q. A Riemannian cubic on G is a curve g(t) that is associated with a solution 
of the HOHP equations (2.4)-(2.6) for L defined in (2.7). If the metric is bi-invariant (that 
is, invariant under both left and right multiplication), then ad^ = — ad, so that the reduced 
Lagrangian becomes u) — | \\u\\^. In this case the HOHP equations (2.4)-(2.6) reduce to 

fi — adt fj,, n = —ii , u = u^, g — g£, £ = u. (2-10) 

These can be rewritten as £ = — ad^ £, together with the reconstruction relation g = g£. 
The solution curves £ are called Lie quadratics in [Noa03]. 

For G = SO (3) we obtain the NHP equation, first derived in [NHP89], 

e=e'xe (2.ii) 

Here we used the identification of the Lie algebra so (3) with M. 3 together with the cross 
product, by means of the isomorphism 



p3 



x)->50(3), Z = {x u x 2 ,x 3 )i->£ 






-x 3 


X-2 







— X\ 


-x 2 


Xl 
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2.3 Hamiltonian structure 

We now demonstrate the symplectic properties of the higher-order dynamics by making 
explicit use of the HOHP variational principle and (2.3). This proof will illuminate a similar 
approach to be used in Section 3 for the corresponding statement in discrete-time mechanics. 
The results can also be viewed in the context of Hamilton's canonical equations on T*(TG). 

The second equation in (2.5) is the constraint relation v — J^, and we are led to define 
the submanifold W c C P&'G that respects this relation. From now on, we assume that the 
reduced Lagrangian £ is hyper-regular, that is, |-§ is non-degenerate. Therefore, the map 
(f) := pr 2 \ w : W c — > T*(TG) is a diffeomorphism, with pr 2 defined in Section 2.1. We remark 
that the Lagrangian £ = | \\ u \\1 f° r Riemannian cubics on manifolds with bi- invariant metrics 
satisfies hyper-regularity since |^| = Id. Here the constraint submanifold W c is determined 
by v = u b . Hyper-regularity is also satisfied when the metric only has one-sided invariance. 

Symplecticity. Assume that equations (2.4) and (2.6) induce flow maps F t : T*(TG) — > 
T*(TG), for all t G [0, a}. These are defined as follows. For an arbitrary x G T*(TG) let 7 be 
the solution to the HOHP principle with initial condition 7(0) = cf)^ 1 {x). Then by definition 

F t (x)=#y(t)). 

Next we show that the flow maps preserve the canonical symplectic form on T*(TG) = 
G x qx 2q*. Let y = (g,£,/i, v) be in T*(TG), and let Vi = {5g u 5^, <5/ij, 5h>i), i = 1,2, be 
arbitrary vectors in T y {T*{TG)). Define r/i = g~ l 5gi. The canonical one-form 6{y) can be 
calculated as 

%)K) = </W> + W£i), (2.12) 
and its exterior derivative, the canonical symplectic form, as 

uj(y){v u v 2 ) = - (5fii,r] 2 ) - + (^2,Vi) + (<^2,5fi) + foi,^]) • (2.13) 

The space of solution curves to the HOHP principle can be identified with T* (TG) by mapping 
a solution 7 to 0(7(0)). Correspondingly, the restriction of Sh p to the space of solution 
curves can be identified with a map s : T*(TG) — > 1R. Namely, if 7 is the solution with 
initial condition 7(0) = _1 (x) for x G T*(TG), then s(x) = 3^(7). Choose an arbitrary 
5x G T X T*(TG), and set ^7(t) = (0 _1 o F t )*(5x) = (Sg^^du^fi^v), the corresponding 
variation through solution curves. We compute from (2.3) and (2.12) 

ds(5x) = [(fi, rj) + (u, 60] a = [{F a ye - 6] (Sx), (2.14) 

with rj as before. The exterior derivative of the generating function relation ds = (F a )*6 — 9 
shows that the symplectic form u is preserved by F a , and therefore by F t for all t G [0, a]. 

Hamilton's canonical equations. We close this section with a brief discussion of Hamil- 
ton's canonical equations, details of which can be found in [CdDll] and [GBHM + 11]. For a 
hyper-regular Lagrangian £ the HOHP equations (2.4) and (2.6) are equivalent to Hamilton's 
canonical equations 

i x uj = dH, (2.15) 
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where X G 3L{T*{TG)) is the Hamiltonian flow vector field and H : T*{TG) ->■ R is the 
Hamiltonian function given by 

H{x) = (//,£) + -^,u(x)), (2.16) 

where |^(£,-u(x)) = z/. As an example, for cubics 

H(x)= 1 -\\vf s + (^0- (2-17) 

As a consequence of (2.15), the Hamiltonian vector field X preserves the canonical symplectic 
form uj on T*(TG). This is an alternative demonstration that F t is symplectic. 



2.4 Momentum maps and Noether's theorem 

In this section we provide a number of necessary definitions concerning group actions and 
momentum maps before discussing the momentum conservation property of the HOHP flow 
map arising from the group symmetry of the Lagrangian. 



General definitions. Consider a left action of a Lie group if on a smooth manifold Q, 

$:HxQ^Q, (h,q)^$ h (q). 



The action, $ : H x Q — )■ Q induces an infinitesimal generator for each ( G f), as the vector 
field (q G 3£(Q) defined by 

$exp H «)(g). 

£ = 

The tangent and cotangent lifts T<3> and T*$ of $ are group actions on tangent and cotangent 
spaces, 



T§ : H x TQ -)■ TQ, (h, v q ) ^ T$ h (v q ), 

T*§:Hx T*Q -> T*Q, (h, a q ) ^ (T$ h _x)*(a ? ). (2.18) 

The cotangent lift momentum map J : T*Q — > f)* associated with the action $ is determined, 
for arbitrary ( G f), by 

(J(-), C> s . x6 = ( • . C Q (q)) T * QxTQ = Oq((t*q(-)), (2.19) 

where 8q is the canonical one form on T*Q and the vector field (t*q G X(T*Q) is the 
infinitesimal generator for the cotangent lifted action (2.18). These cotangent lift momentum 
maps are an important special case of the general concept of a momentum map of a group 
acting on a Poisson manifold [MR03]. 



Burnett, Holm, Meier Higher-order mechanics on Lie groups 



11 



Noether's theorem. We proceed with a version of Noether's theorem for systems with 
symmetry. The action $ of a Lie group H on G induces an action on the space of curves 
C(G) and thus on both TG and T^G. Let us consider a Lagrangian L : T&G ->• R, with 
the symmetry 

LWS), (2-20) 



L ([$, 



l(2) 



for all h G H and any [g]^ G T^G. For consistency we continue to assume the Lagrangian 



is left-invariant, however the assumption is not essential for this section. For a curve g(t) G G 
consider the variation induced by the action of H, 5^g = Cg{q) for a £ G P). Then 

5 ( [ a L([g])=0. (2.21) 
Jo 

If g(t) is a solution to Hamilton's principle (2.1), then by Theorem 2.1 there exists a lift 
7(t) G P {2) G satisfying the HOHP equations (2.4)-(2.6). By (2.14) and (2.21), 

= ds(6 ( x) = 6(6 c x(a)) - 6{6 c x(0)) = 6(( T *(T G Ma))) - e(( T * {TG) (x(0))), 

where x(t) = <f>('j(t)) G T*(TG). Using (2.19) we identify the cotangent lift momentum map 
J : T*(TG) — > \)* associated with the action T$ : H x TG — > TG and find the conservation 
law 

J(x(t)) = J(x(0)) for all t G [0, o]. 



Example. Let H = G and consider the left action L of G on itself, 

L : G x G — ^ G, L h (g) = hg. 

This is by definition a symmetry for any left-invariant Lagrangian. The tangent lift TL is 
given by 

TL : G x TG, TL h (g,£) = (hg,0, 

where as usual we used left trivialisation, TG = G x q. The cotangent lift momentum map 
associated with TL is 

J:T*(TG)^ *, (g,Z,ti,v)^Ad* g - lf i. 

Therefore, J is conserved along solutions of the HOHP equations (2.4)-(2.6). 

The Hamiltonian approach to Noether's theorem is discussed in [MR03]. It is however 
the variational viewpoint that will expedite similar results for the variational integrators we 
develop next. 



3 Discrete-time mechanics 

Continuing with the Hamilton-Pontryagin approach of the previous section, we now address 
the issue of geometric discretisation. We proceed systematically with a discrete variational 
principle and, by considering end point terms, we establish the symplectic and momentum- 
preserving properties of the resulting numerical scheme. 
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3.1 Runge-Kutta Munthe-Kaas (RKMK) methods 

In this section we obtain a discretisation of the reconstruction equations (2.6). This will be 
achieved by combining a standard Runge-Kutta (RK) method to integrate £ = u with a 
Runge-Kutta-Munthe-Kaas (RKMK) method to integrate g = g£. We shall introduce both 
types of methods in sufficient detail for our purposes here and refer to [HLW06, IMKNZOO] 
for more complete discussions. 

Runge-Kutta (RK) method. An s-stage RK method with step size h may be employed 
to numerically integrate first-order ordinary differential equations of the form y = f(y,t) 
with y(t ) = y . The task is to find solutions Y l for i — 1, . . . s that satisfy 

s 

Y^yo + hJ^^jfiY^to + Cjh). 
i=i 

The numerical estimate of y(to + h) is then given by 

s 

yi=yo + hY^ bif(Y\ t + ah). 
i=i 

The Y % are known as internal stage variables and the real coefficients can be presented in a 
Butcher tableau, 



Cl 


an ■ 


■ a u 


C s 


a sl . 






bi . 


■ b s 



These coefficients must obey q = Ylj a ij an d Yli h — ^ f° r a n RK method of at least first 
order accuracy. A discussion of order conditions can be found in [HLW06]. 

An s-stage RK method for the reconstruction equation £ = f(£,t) = u(t) starting from 

f (*o) = £o is thus 

s s 

& = £ + h-Y i aijV j , Zi = Zo + hY,biV\ (3.1) 

j=i i=i 

where := u(t + Cjh) and & is our estimate of £(t + h). 

Runge-Kutta— Munthe-Kaas (RKMK) method. We now turn to the reconstruction 
equation g = g£. In the RKMK method, one begins by rewriting the reconstruction equation 
on q, to which one then applies an RK method. Consider an approximation to the Lie 
exponential r : q — > G. Let r be an analytic diffeomorphism in a neighbourhood of with 
r(0) = e and t(—T])t(t]) = e for all r] G Q. We define the right-trivialised differential dr^ at 
£ e S as, 

dr e : ->■ 0, V ^ (^t(77))t(0 _1 , 
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with inverse dr^ . We record the formula 

dT i \ V ) = dr:^Ad T{ ^ )V ), (3.2) 

which can be obtained by differentiating t(— £)r(£) = e. For small times t a solution of 
g = g£ with g(0) = go can be parameterized as g(t) = gor(Q(t)) for a © G C(g). Then 

9Z = 9oT(Q)£ = g = g (dT e (Q))T(Q), 

and therefore, using (3.2), 

Q = f(Q,t) = drZ 1 m- 
An RK method for this equation, with starting point ©(to) = 0, is 

s s 

0* = ]T flif^, 6! = h hdrZl ei ~\ (3.3) 
j=i i=\ 

where S J := £(io + Cjh). Define also G l = goT(hQ l ). Equations (3.3) together with 

9i = 9or(&i) (3.4) 

constitute an s-stage RKMK method with step size h. Equation (3.4) is our estimate for 
g(to + h). The following theorem can be found in [HLW06]. 

Theorem 3.1. Assume the underlying RK method in RKMK has order p and r is a p th - 
order approximation to the Lie exponential. If one truncates dr' 1 at order p — 2, the RKMK 
method is still of order p. 

For p = 2, therefore, equations (3.3) can be simplified by approximating dr~ x by the 
identity map. The resulting RKMK method is of second order, 

s s 

T- 1 (g Q - 1 G i ) = hJ2a ij E\ r{g, l gi ) = h^b0 . (3.5) 

j=l i=l 

Discrete reconstruction equations. Our final method for effectively solving (g~ l g)' = u 
comes in the two stage form achieved by combining the above methods, vector space RK and 
RKMK. In (3.5), 5* are exactly ^(to + c «); i n (3-1) the S l are thought of as being approximately 
£(to + Ci). The combined method feeds the calculated in (3.1) into (3.5). We now proceed 
by inserting the discrete reconstruction equations into a discrete HOHP principle. 

3.2 Higher-order Hamilton Pontryagin Integrators 

Discrete path space. Our goal is to approximate the reduced higher-order Hamilton- 
Pontryagin principle (HOHP) in (2.2). This is achieved by replacing the integral with a sum 
over discrete time points t%. = kh for k — 0, . . . , N, where h is the time step. The continuous 
variables are replaced by the s-stage RKMK variables discussed above and the kinematic 
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constraints by the RKMK equations. A tentative choice of discrete path space is therefore 
given by the set of maps 

{7:{t k }to^T*(TG)x(pWG) S } 
which, by using the variables 

l(h) = (g k , to, "k, {Gt S* fc , V£, /4, vl\ N s=l 

can be identified with the set T*(TG) N+l x (P^G) S ' (N+1) . The discrete cost functional that 
we will define below in (3.6) will not depend explicitly on Hq 1 uq or G^y, Sjv, Vn, I^n, ^N- 

We 

therefore choose to omit these elements and define the discrete path space Cd as follows, 

C d := TG x (T*(TG)f x (P^G) N - S 

= (G x q) x (G x q x 2 5 *) N x(Gx2gx 2q*) n - s . 

We shall use the following index notation to denote elements, 

(go, £o)> (9k, €k, Afe: Vk)k=i , ( G l, s l) Yk> f4>, u k)*'k=i o 



Discrete variational principle. We now give the discrete HOHP integral. Define the 
functional S : Cd — > R as 



N-l 



k=0 



s 1 s 

1=1 

1 s 



f(/i fe+ i,-r ^V+i) - ^6iHJ. 

i=l 



(3.6) 



Let 7 G Crf and £7 G T 7 Cd. Also denote % = g k l 5gk and r^jL = (G l k ) 1 SG t k . Moreover, set 
^fc+i = Sj=i ^j^j an d ©fc = Sj=i a u'"i- Furthermore, it will be useful to define 



u Q = v x 



as well as 



s s 
t=l i=l 



A** = ( dT -hr. k T ^k, k = l,...,N. 



>0 



(3.7) 
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In the following equations, the index k takes range 0, ... N — 1 unless explicitly stated oth- 
erwise. Discrete variations of 5* now give, 

SS = h J2 ( [J2 (\( dT -h0])*ri°> Vk ) + { hi Jt + \ Vk ~ 2 a ^A - 6 i#*+i, 5Hj 



<5w 

3 

N-l 



k=l i i 

+ (jjln, Vn) + (vn, $£,n) - (/io, Vo) - (vq, 8 to) ■ (3-8) 

A curve 7 G Cd is a solution of the HOHP principle, if 5S = dS^^) = for all vectors £7 
with fixed boundary points g , g^ and velocities £ 0) £iv, that is r] = t]n = <5£ = 8£n = 0. We 
set the notation = J|(Hjj., V fc 4 ) and obtain the RKMK reconstruction equations 



(3.9) 



i=i j=i 

6+1 = £* + X] ^'^fc ' 9k+l = 9kT ( h Y b i E k 
5=1 J=l 

the auxiliary equations 

/4 = , v k +i = v k - Y "*> ( 3 - 10 ) 

j 

the discrete equations for the Ostrogradsky momenta, 

§ = « + , - ±4 , | = £ (f - l) 4 - ^ (3.11) 

and the discrete version of the Euler-Poincare equation, 

^fc+i = (^ Hfe+1 )*(^H fc+1 )>fe- (3.12) 



3.2.1 Geometric properties 

Symplecticity. A solution 7 e is said to have initial conditions (go, £0, /^O; ^0) £ T*(TG), 
upon using the definitions of and /^o given in (3.7). Assume the Lagrangian t allows the 
internal variables G\, E, l k , V k l , fi k , v\ to be eliminated from the equations, which subsequently 
induce a discrete flow map F : T*(TG) — > T*(TG). That is, for any solution 7 the flow map 
satisfies F k (go, £0, fJ>o, ^0) — (g k ,£,k, ^ k ,Vk)- We express the discrete HOHP cost functional 
in terms of initial conditions, as follows: Let s : T*(TG) — > R be defined by s(z) = Sl'j), 
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where 7 is the solution curve with initial conditions z £ T*(TG). A variation 5z G T Z T*(TG) 
corresponds to a variation £7 G T^Cd, and we compute the generating-function relation, 

ds(Sz) = (li n ,Vn) + (vn,S£n) ~ (jio,m) ~ (^0,^0) = [(F*f0 - 0] (Sz) . (3.13) 

Upon taking the exterior derivative, one finds that the symplectic form u is preserved by 
F . In particular, setting N = 1 shows that the discrete flow map F is symplectic. 

Momentum preservation. The argument continues along the same lines as in Section 2.4. 
Let us assume we have a (left-invariant) Lagrangian L : T&'G — > M. and an action $ of a Lie 
group H on G. Note (G\, HJ., V fc l ) G G x 2q = T^G. For the group action to be a symmetry 
of the discrete dynamics we require that the Lagrangian is group invariant, as in (2.20), and 
the group induced variations on (6^,5^, V£) G T^'G preserve the reconstruction relations 
(3.9). That is, for ( G f) and h(e) := exp if (eO the induced variation (G\(e),E\(e),V£(e)) 
still allow (3.9) to hold for some {^(e), 6fc( e )}fc=o- Similar to the continuous case, if all the 
underlying RK relations between the variables in Cd are implicitly assumed, then a solution 
(G£,Hjj., VJ?) to Hamilton's principle 

5j>£(HLV;> = 0, (3.14) 

k,i 

can be lifted to a solution 7 of the discrete HOHP principle. Note that here variations satisfy 
fixed end point conditions 5go = 5g^ = and <5£o = 5£jv — 0. The group symmetry implies 

* f 5>*(5iX)* = 0- (3.15) 

For a variation ^7 respecting (3.9) this implies S^S = 0, for S as in (3.6). From (3.13) we 
can then conclude that the cotangent-lift momentum map J : T*(TG) — » f)* associated with 
the action T$ : x TG — > TG is preserved by the HOHP integration scheme, 

= ds(5 c z) = 9(Ct*(tg)(zn)) - 0(Ct*(tg)(zo)) = J(^v) - J(«b), 
where we wrote z G T*(TG) for the initial condition of 7 and z^ := F N (z ). 

Example. As in the continous-time case, let H = G and consider the cotangent-lift mo- 
mentum map associated with TL : G x TG — > TG, the tangent lift of the left action of G on 
itself. We recall that 

J:T*(TG)^Q*, {g,Z,»,u)^Ad* g - lf i. 

By the above discussion, J is conserved along solutions of the discrete HOHP equations 
(3.9)-(3.12). That is, := Ad*-i is conserved, J^+i = J&. 

We remark that this conservation law can alternatively be obtained directly from the 
HOHP equations and a property therein which relates them to the Euler-Poincare equations. 
It follows from (3.2) that for any rj, £ G Q and //£$*, 

((dr^Y (n), V ) = {li,drZl{r))) = (fi, dr^(Ad T(0 v )) = (Ad^drf ^(j*)^) • 
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Therefore 



Ad; 



(3.16) 



For a solution of the discrete HOHP equations, (3.12) implies 



(J>k+i 



(drl 
Ad! 



-i 

hS k+i 



)*(dT hSh+1 )*n k = Ad* T{hEk+i) y, k 




■r(hSk+i) 



Hence, 



Ad* 



i Mfe+i 



Ad* 



'Sfc+i 



which is the conservation law for the discrete momentum, Jfc+i = J*,. 
3.2.2 Implementation 

Our primary application of the HOHP scheme addresses the Riemannian cubics on 50(3) 
with a bi-invariant metric. We implement initial value problem solvers using the implicit 
Euler method and the Stormer-Verlet method. 



Fig. 3.1: An illustration of the Stormer-Verlet method applied to the initial value problem of Riemannian cu- 
bics on the rotation group, presented in (3.17)-(3.19). In this figure, a given element of the group corresponds 
to a point on the radial line along the rotation axis. The distance of the point from the origin represents 
the rotation angle measured in a right-handed coordinate system. The center and the boundary of the outer 
sphere therefore both represent the identity matrix. The displayed curve results from numerically integrat- 
ing system (3.17)-(3.19) over a period of 2ir. The chosen initial conditions were go — Id, £o = (—6,1,0), 
vq = (0, 0, 6) and \i§ = (0, 36, 0), which belong to certain 2-7r-periodic solutions found in [NoaOl]. 
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Numerical schemes. We consider two instances of the general RKMK scheme, an implicit 
Euler method and a Stormer-Verlet type method given respectively by, 



1 1 ° 
— and 1 





1/2 1/2 



1/2 1/2 

For the Riemannian cubic Lagrangian, £(£,u) = ^\\u\\ 2 , the Euler method gives 

/i fc+ i = {drZl^ k+x T{dr Kk+1 )*n k , 
v k +i = v k - h(dT-h£ k+1 )* A*fe+i , 
9k+i = 9kT(h£ k+1 ) , = £ k + hv\ . 

This method is explicit and straight-forward to implement. For Stormer-Verlet the result is 
similar but implicit 

= (dTZl Sk+i )*(dT h3k+1 )*n k , E k+1 = i (f fc+1 + £ k ) , (3.17) 
Vk+\ = v k - h(dr hSk+1 )* fi k , (3.18) 

gk+i = g k r(hE k+1 ) , £ k+ x = £ k + h \v k - ^(dT hSk+1 )* /ikj ■ (3.19) 

The most notable differences are the appearance of the average velocity and the second 
order approximation in the equation. To implement this method we can apply a fixed 
point method to find with the function 



fk{0 =£-tk-h(v h + 7i (dTh^^ fi k ^j 



An example of a map r is the Cayley map, cay : so(3) — > SO(3), given by 

cay(X) = U-X/2\ 1 [e + X/2 s j, d cay^ Y = (e - X/2\ 1 Y (e + X/2\ 1 (3.20) 

where e is the identity matrix. With the Cayley map we can implement the Stormer-Verlet 
method, say, on 50(3) to produce plots such as in Figure 3.1. Note that both algorithms 
give discrete flow maps F : T*TSO(3) — > T*TSO(3) which are symplectic and momentum 
preserving, as discussed in Section 3.2. 



4 Application to trajectory planning 

In this section we discuss an application of the discrete HOHP variational principle to a 
second order trajectory-planning problem, where the interpolating curve evolves by means of 
a group action. One aims at finding a curve that passes near a series of given target points at 
prescribed times. A precise definition of the problem will be given below. An earlier inves- 
tigation of the same problem in [GBHM + 11] found that optimal curves satisfy second order 
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Euler-Poincare equations between time nodes. At the time nodes the (otherwise conserved) 
momentum experiences a kick that is related in a simple fashion with the optimal mismatch 
between interpolating curve and target at the respective node, leading to a piecewise con- 
stant momentum. The goal of the present section is to show that this geometric characteristic 
carries over to a discrete-time algorithm. 

4.1 Continuous-time trajectory planning 

Background. Our motivation to study this type of trajectory-planning problem lies with 
potential applications in computational anatomy, the modeling and quantifying of diffeomor- 
phic shape evolution [MTY02, MY01]. Usually one seeks a geodesic path on the space of 
shapes between given initial and final data. This approach can be adapted for longitudinal 
data interpolation, interpolation through a sequence of data points, by piecewise geodesic 
curves. A class of higher-order generalizations providing smoother solution curves than the 
piecewise geodesic ones was studied in [GBHM + 11] in the finite dimensional setting. 

Here we first reformulate the problem in terms of the HOHP principle, which is amenable 
to the discretisation introduced in Section 3.2. 

Problem formulation. Consider a left representation of the Lie group G on a vector space 
V with norm ||.||y, 

GxV^V, (g,I)^gL 

The problem under consideration is defined as follows: Given a Lagrangian i : 2q — > R, 
a, tx, ■ ■ ■ , ti G M ; To, It x , . . . , It x G V , and £q G q, minimise the functional 

S := T £(£, u) + (fi, g~ l g -^(i/.j-^+l^ \\g'\U)T G - I u f y (4.1) 
Jq 1 i=i 

subject to the conditions £(0) = £q an< ^ fi'(O) = e - Variations are taken amongst curves 
(g,£,u, fjL,p)(t) G P( 2 >G, where we require that g G C^QO, 1]) and g\[t u t l+1 ] G C°°([ti,ti+i]). 

This type of trajectory-planning problem is familiar, for example, from computational 
anatomy, where one would typically think of To as a template that is deformed by a curve of 
diffeomorphisms g -1 (t), in turn generated by the time-dependent vector field £(£). At times 
ti the curve passes near the given targets I t ., the parameter a determining the proximity 
of the passage. In that case the Lie group G is infinite dimensional. It is clear that this 
type of trajectory-planning problem is also relevant in numerous finite-dimensional situations 
whenever a path g(t)To generated by transformations g(t) is required to pass by prescribed 
target points I t . at given times tj. Here we focus on the latter case. 

Euler-Lagrange equations. We assume that the norm on V is induced by an inner prod- 
uct (•, • ) v and define the isomorphism b : V — > V* by (u, v) v = (w b ,t>). We denote the 
cotangent-lift momentum map associated with the action of G on V by o : V x V* — > g*, 
that is (see (2.19)), 

(/ou;,O fl . XB = <w,£vOO)y.xv> for all / G V, u G V*, £ G Q. (4.2) 



Burnett, Holm, Meier Higher-order mechanics on Lie groups 



20 



We compute the Euler-Lagrange equations governing any minimiser of S in analogy with 
earlier calculations done in Section 2.2. The only complication arises from taking variations 
of the penalty term in (4.1). Setting rj := g~ 1 5g ) we obtain 



/ 



5 i E H^r'To - hf v = - A E ( (ai^-'To - 4)' . (v(U)v) (g(U)T 



i=l i=l 



E ((^r lT o) o (^(^)- i t - ij b ,^ 



i=l 



Taking variations of (4.1) and requiring 5S = yields the HOHP equations (2.4)-(2.6) on 
open intervals We recall these equations as 

fx = ad ? n, n= — -v, v = — , g = g£, £ = u. (4.3) 

Moreover, writing tf for the limits from above and below respectively, 

"(tt) - "fa) = 0, (4.4) 

- M*D + ^ fafe)" 1 ^) o {giur'n - I u f = 0, (4.5) 



for i = 1, . . . I — 1, and 

K*0 = 0, (4.6) 

Kti) - l~ 2 (g{tiT l T Q ) o (g{tiT l T Q - i M y = o. (4.7) 

Momentum kicks. Optimal curves satisfy the HOHP equations (4.3) between time nodes. 
While the Ostrogradsky momentum v is continuous in time (4.4), the Ostrogradsky momen- 
tum fj, experiences a discontinuity (kick) at the time nodes. The kick is related to the optimal 
mismatch between interpolating curve and target points at the respective node, (4.5). The 
final values at t\ of v and \x are given in (4.6) and (4.7). The prescription of target points 
I t . breaks the symmetry of the problem, which is reflected in the fact that the momentum 
J := Ad i/i is no longer preserved. Indeed, J is piecewise constant, with discontinuities at 
the time nodes given by 

J(tt) - J(tr) = —Ad^ ((g&r'To) o {git^T, - I ti y) . (4.8) 



4.2 Discrete-time trajectory planning 

As in Section 3.2 we replace the integral in (4.1) with a sum over time points % = kh for 
k — 0, . . . , N, where h is the time step. Again the continuous variables are replaced by the 
s-stage RKMK variables and the kinematic constraints by the RKMK equations, see Section 



Burnett, Holm, Meier Higher-order mechanics on Lie groups 



21 



3.1. Let the target times be U = Nih for i = 1, . . . I with iVj = N. For convenience we also 
define Nq := 0. The discretisation of (4.1) is given by 



N-l 



k=0 



s 1 s 

i=i " i=i 



+ (/i fe+1 ,-r 1 (( 7fc 1 ( 7fc+1 )-^^) (4.9) 



i=l 



1 S 1 

1 i=l ' »=1 



|2 



Discrete Euler-Lagrange equations. For the variations of the penalty term we obtain 
1 1 1 1 

6 ^ E ll^ T ° - 7 <* II v = -^2 E ( W T °) ° te T ° - J b ' ^ 



i=l i=l 

where we set rjk '■= g^^dk- As in (3.8) we define 

= (drZl Sk Yfi k , k = 1, • • • , N. 
Taking variations of (4.9) and requiring SS = gives, for k = 0, . . . N — 1, 



(4.10) 



s s 

2* = a + E ' G fe = ^ fcr E a ^ 

3=1 J=l 

s s 

&+i = 6; + /i E & ^ ' ^ fc+i = ^ r E ^ E 

i=i j=i 



(4.11) 



i4 = > 

^fe+i - ^ + E ^ = 



A; / j \* ^ i A: 

— = (rfr_, Sfc+1 ) ^ ; 



E 



Moreover, for interior indices k ^ Ni (i = 0, . . . /), 

J/fc - (^H fe+1 )*( rfr -^ fc+1 )>fe+l = 0, 

and for node indices N{ (i = 1, . . . I — 1), 



(4.12) 



^ - («*r^. +1 )*(dr_ fcaWi+1 )>^ + i - ^ (^T ) o (gj+T - I k )" = 0. (4.13) 



Finally, 



v N = 0, 



(7 



(4.14) 
(4.15) 
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Discrete momentum kicks. Note that equations (4.11) are a subset of the HOHP equa- 
tions of Section 3.2. The last equation in (4.11) is the discrete version of (4.4). The discrete 
Euler-Poincare equation (4.12) completes the discrete HOHP equations for interior indices. 
This corresponds to the continuous-time equations (4.3), which hold on open intervals. At 
time nodes iVj, the update formula for fi^+i acquires an extra term, described by (4.13). 
This is analogous to the discontinuity of fi(t) seen in (4.5). The exact correspondence be- 
tween equations (4.14) and (4.15) with (4.6) and (4.7) is obvious. The kicks of the discrete 
momentum := Ad*_i fik are computed from (4.12) using (3.16), 

J Ni+1 - J Ni = ~ AdJ-j ((^jTo) O (gj*T - I ti ) b ) , (4.16) 

for i — 1, . . . I — 1. This is the exact analogue of (4.8). If k is an interior index, then J^+i = 
This discussion shows that the geometric behaviour of continuous-time momentum carries 
over to the discrete-time setting described above. 

4.3 Implementation via shooting method 

In this section we address the practical problem of minimizing the discrete functional (4.9). 
Since any minimizing curve is known to satisfy the shooting equations (4.11)-(4.13), the 
search can be restricted to the space of solutions to these equations which results in a problem 
of reduced dimensionality (see also [VRRC11]). The functional (4.9) is accordingly written 
as a function of initial momenta /io and uq, J : 2q* — > R. A gradient descent algorithm is 
then employed to carry out the minimization of J . We discuss in some detail an efficient 
way of computing the gradient Vj7 and present numerical simulations for the case of SO (3) 
acting on M 3 . 

Shooting method. Explicitly, J : 2g* — > R is 

J(jM>, Uq) = T u)dt + lk -1 (^) T o - 4 II v > ( 4 - 17 ) 

^° 1=1 

where equations (4.3)-(4.5) are implied, with initial conditions g(0) = e, £(0) = £° (recall 
that these are prescribed by the interpolation problem), /i(0) = /io and z/(0) = u . That is to 
compute j7"(/i(b ^o); integrate (4.3)-(4.5) with the given initial conditions, then evaluate the 
right hand side of (4.17). 

In discrete time one proceeds in analogous fashion, defining : 2q* — > R as 

AT-l s i 

Mho, "o) = hJ2J2 b ^ ^) + 2^2 E \\9nIT - h\\ 2 v , (4.18) 

k=0 i=l i=l 

the discrete shooting equations (4.1 1)— (4.13) being implied. We then employ a gradient 
descent method to minimise Jd- By construction, the corresponding solution curve obeys 
momentum map conservation for interior indices and has momentum kicks (4.16) at node 
indices. 
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4.3.1 Computation of the gradient 

We now discuss the computation of the gradient VJd via adjoint equations, the concept 
of which we first expound in continuous time before presenting the discrete case. For this 
discussion we choose the bi-invariant cubic Lagrangian £(£,u) = ^\\u\\ 2 . 

Continous time. For simplicity, let / = 1. Introduce adjoint variables P°, P 1 , V°, V 1 , V 2 G 
20* x 3q =: A used to explicitly enforce the implied constraints in the definition of J . Take 
S : C(P^G xA)^lto be 

5(a) = [ 1 ^\\u\\ 2 s dt + ^Wg-'it^To - I, 1,2 



o 



illy 



+ J 1 <p°, g- l g - £) + (p\ £ - + + v°) (419) 

+ (ft - ad* fi, V 1 ) + (v - u\ V 2 ) dt . 

Note that if a, the projection of the curve a G C(T*(TG) x A) onto T*(TG), satisfies 
the HOHP cubic equations (2.10) with initial condition 5(0) = (go, £o> A*o> ^o) then S(a) = 
J~(fio, Vq). Taking variations while keeping (g(0), £(0)) fixed (the prescribed part of the initial 
conditions) we find that if such a curve a obeys terminal conditions, 

P°(h) = l~ 2 (gihrn) o (gihm - i tl ) b , (4 2Q) 

P 1 (t 1 ) = 0, V°(t 1 ) = 0, v\t 1 ) = o, 

and adjoint equations 

P° = ad*P°, Pi = -P° + ad yl ^, 
V^ = u -(P l f 1 V 1 = V° - ad 5 V 1 , 

then the gradient VJ is given by 

V W J = -V^O), V U0 J = -V°(0). (4.22) 

In practice, one integrates (2.10) forward and (4.21) backward in time, with initialisation at 
t = tt given by (4.20). 



Discrete time. For the purposes of implementation we are primarily interested in the 
discrete case and we now derive adjoint equations consistent with our discretisation of the 
forward shooting equations. We now allow for I > 1, but for simplicity restrict to the Euler 
type HOHP integration scheme for which the equations with momentum kicks (4.11)— (4.13) 
can be written as, for k — 0, . . . , N — 1, 

9k+i = 9k r {h^k+i) > ^k+i = £fc + hut , ik+i = ik + hut , Uk — v k , ,^ ^ 

Vk+i = Vk~ hfik+i , {drj£ Yfrk+i = {drZ^Tfik + ®k{9k) , 
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where <& k (g) '■= for k 7^ Ni, i = 1, . . . , I and 

<$> k {g):=-- 2 {g- l T )o(g- 1 T -I t )\ for k = N i} i = l,...,l. (4.24) 
o~ z 



The discrete version of (4.19) is 



N - l h, 



Sd(a d ) = ^2^\\u k \\ g + (P k+1 ,r l (g k 1 g k+ i) - hE k+1 ) + (.F& f i>&+i - & - hu k ) 

k=0 

+ (P k+1 , - £ fe - hu k ) + (z/fc+i - z/ fc + hfa+i, V k ° +1 ) 



(4.25) 



1 ' 2 



i=i 



Notice that the discrete adjoint space Ad has an extra variable P 2 6 g*. The resulting discrete 
adjoint equations consist of the terminal conditions 



h^tf) 1 N - N\yN ) 1 

Pi = -DtVl = , V# = , vi = 



TV u ' ' AT 

and the discrete adjoint equations 

P° k = (dr^ k )* ((dT^ k+i yP° k+1 - A k (g k , V k l +1 ] 

Pi = hP° k+l + Pl +1 + D'V^ - D k V k , 
V k ° = V k ° +1 + h{Pl +l + hP° k+1 - u k ) , 
V, 1 = dr h ^ k (drZ^V^ - hV k °) . 

In the above we used the abbreviation := where the maps Df^ '■ — > Q* are 

defined by the relation 



d 

de 



( dT ±h(Z+ep)) ^ a ) . = ( D t» a >P) S *x 3 ' fOT a11 P G S • ( 4 - 27 ) 



Moreover we defined A k : G x $j* — > q* as 

(A(<7,^),^) := ($fc(<7),?7> + ^ 



($ fc (^exp (erj)),V) 



e=0 



for all 77 G 0. Then the gradient of the functional defined in (4.18) is given by 

V^J d = ~V \ V U0 Jd = -V °. (4.28) 
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4.3.2 Interpolation on S 2 

We present simulations for G = SO (3) and V = M 3 , where the group action is by matrix 
multiplication of vectors. In this setting the o operator of (4.2) becomes the standard cross 
product for I,u G M 3 , 

/ow = /xwel 3 =so(3). 

For our purposes we take the starting point To = (1, 0, 0), initial velocity = (0, 0, 1), the 
node times t{ = |z for i = 1, . . . , 5 and target points 




Note that all target points were chosen to have unit length so that we effectively perform 
interpolation on S 2 . We refer to Figure 4.1 for an example interpolant and Figure 4.2 for an 
illustration of the momentum behaviour. 



Fig. 4.1: Discrete-time trajectory planning on the 
sphere. The curve shown is a minimiser of the cost 
functional (4.9) obtained numerically by following the 
procedures outlined in this section. Namely, the cost 
functional was restricted to solutions of the shooting 
equations (4.11)-(4.13), thus reducing the problem to 
the space of initial momenta, where a gradient descent 
was performed. The initial and target points were cho- 
sen as given above with tolerance parameter a — 0.025. 
The colours represent the local speed along the curve 
in SO(3), that is ||^|| (red is large, white is small). 



5 Summary and Outlook 

This paper has discussed a structure-preserving numerical integration scheme for a class of 
higher-order mechanical systems. The structure-preserving properties arose as a a conse- 
quence of the variational nature of the discretisation method. 

Our starting point was a continuous time higher-order Hamilton-Pontryagin variational 
principle, where second-order kinematic constraints were enforced via Lagrange multipli- 
ers. These kinematic constraints were transferred to the discrete-time setting by combining 
Runge-Kutta and Runge-Kutta-Munthe-Kaas methods, in a similar way to the first-order 
treatment in [BRM09]. The resulting discrete flow maps were shown to preserve the symplec- 
tic form of higher-order mechanics and to respect momentum conservation in the presence 
of group symmetries. We proceeded with an application to a trajectory-planning problem 
familiar from inexact template matching in computational anatomy. The continous-time 
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Fig. 4.2: Momentum norms. For the interpolating discrete cubic of Figure 4.1, the plot shows the norms of 
the momenta /i^ and v^. The norm of [i^ displays momentum kicks at node indices and exact conservation by 
discrete coadjoint motion for interior indices, as in equations (4.12) and (4.13). The norm of V}. demonstrates 
continuity as found in the last equation of (4.11). Both graphs respect terminal conditions (4.14) and (4.15). 

problem was presented in the higher-order Hamilton-Pontryagin framework, which led to a 
particularly clear exposition of the geometric properties of solution curves. These proper- 
ties were finally shown to be inherited by the higher-order Hamilton-Pontryagin integration 
scheme. 

The treatment in the present paper was focused on mechanics based on Lagrangians 
depending on velocity and acceleration. It is however clear that a generalization to third, 
and higher, order follows in a straightforward manner. Other directions for further research 
include a rigorous analysis of the level of accuracy of the higher-order Hamilton-Pontryagin 
integration schemes. A related goal would be to include Runge-Kutta-Munthe-Kaas methods 
with a higher order of accuracy by choosing a larger truncation index in the sense of Theorem 
3.1. 

The type of trajectory-planning problem treated here is relevant in a number of finite- 
dimensional situations, whenever a path g(t)T generated by Lie group transformations g(t) 
is required to pass near prescribed target points I t . at given times £j. We have treated the 
case of rotations, in which the Lie group is 5*0(3). However, one may envision applications to 
many other problems involving different Lie groups. For example, one could imagine applying 
these methods in designing particle-beam optics, in which the rotations would be replaced by 
symplectic transformations acting on phase-space moments of the beam distribution function, 
such as its emittance [Drall]. 
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